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Abstract. We derive the Bosonic Dynamical Mean-Field equations for bosonic atoms in 
optical lattices with arbitrary lattice geometry. The equations are presented as a systematic 
expansion in l/z, z being the number of lattice neighbours. Hence the theory is applicable 
in sufficiently high dimensional lattices. We apply the method to a two-component mix- 
ture, for which a rich phase diagram with spin-order is revealed. 
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£Sj 1. Introduction 

-f-, Bosonic atoms in optical lattices, described by the Bose-Hubbard Hamiltonian HI 12, 

^vq can easily be brought into the strongly interacting regime by increasing the intensity of the 

C\) laser beams forming the optical lattice [3 1. At commensurate fillings this is reflected in the 

I/"} formation of a Mott Insulating state. This transition has been detected experimentally [4|, 

£ — as one of the first quantum phase transitions in the strongly correlated regime in the field 

^^ of cold atoms. 

In the case of a mixture of different bosonic components, either of different hyperfine 
states of the same atom or different species, additional spin order can exist in the Mott 
k* insulating state, thus further enriching the phase diagram J5] |6l [7] |8l [9] [10l QT| [12) . 

Whereas to describe the single-species Mott insulator transition it is sufficient to use 
Gutzwiller Mean-Field Theory lfT3lfT4l . which uses the local superfluid order parameter 
as the mean field, to capture spin order one needs more sophisticated methods. The reason 
is that the superfluid order parameter vanishes in the insulating states, effectively decou- 
pling the sites, such that all possibilities for spin-order are degenerate. Here we formulate 
Bosonic Dynamical Mean-Field Theory (BDMFT) (BJITTlEDiniEDEDIlD), which is 
able to describe long-range spin order, because in addition to the Gutzwiller Mean-Field 
term describing the superfluid-insulator transition, also a dynamical coupling is present, es- 
sential for the description of virtual hopping processes in the insulating states. The theory 
is presented here as a systematic expansion in 1 /z, z being the lattice coordination number. 
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A generic bosonic mixture consisting of N flavors in a sufficiently deep optical lattice 
is described by the following Bose-Hubbard Hamiltonian 

(1) "H - - ^ J a (b\ a b ja + h.c.) - ^ n a h ia + - ^ U ay h ia (n iy - 6 ay ) 

{if), a i,a i,a,y 

Here Greek indices (a, y) denote the flavors and Roman indices (i, j) denote the lattice 
sites, such that b denotes the annihilation (creation) operator of a boson with flavor a at 
site i. We have introduced the number operator h la - b. b. . The tunneling parameters 
J a can be different for different species, which is relevant when the atoms have a different 
mass or a different coupling to the laser beams. Also the chemical potentials fi a can be 
different, in order to allow independent tuning of the densities of the different components. 
The inter- and intraspecies interactions are encoded in U ay . 

2. Methodology 

In deriving the BDMFT equations we consider the limit of a high- (but finite-) dimen- 
sional optical lattice lfTTll2TI . Hence z, the number of neighbouring lattice sites, is large 
and 1 lz can serve as an expansion parameter. In order for the total energy per site to con- 
verge in the limit z — > °o, the tunnelling amplitudes then need to be rescaled as J a = J' a /z. 
This introduces the small control parameter 1 /z into the Hamiltonian. At the final step in 
the derivation we will return from J' a to J a . We thus present the BDMFT equations fol- 
lowing a systematic expansion in the small parameter \jz ifPTl . This is different from the 
original proposal ITT31 [181 [191 . as well as from a recent derivation 111711201 . although the 
final equations in all cases (for finite dimensions) coincide. 

2.1. Effective Action. Using the conventional coherent state representation (see e.g. (22)), 
the partition function Z is expressed as 



/' 



(2) Z= D[b*]D[b]exp(-S[b*,b]/fi) 



where we use the shorthand notation jD[b*]D[b] = f Y\i, a D[b* a ]D[b ia ], and the b ia {T) 
are complex- valued fields PI The action S [b*, b] is given by 



(3) 



5 [b*,b] = f dr\ J] b] a {r) L^- A Mt) - 2 7 [K(T)b ja (T) + ex.] 

{ ',a ' ' (ij),a 

+ \Yu U ay b] a (T)b ia (T) [b* y (T)b iy (T) - 6 oy ] 1 . 
i,a,y j 

Following the same 'cavity' method used to derive the fermionic DMFT equations ED . 
we now consider a specific site and call it site 0. We split the action into three parts: Sq 
contains those terms exclusively related to site 0, S (0) contains those terms not incorporat- 
ing site 0, and AS contains the terms connecting site to the other sites. The system with 
site excluded is called the cavity system. AS is thus given by 

(4) AS [b* , b] = - i J drj] J' a [b;jr)bja(T) + ex.] , 



(0j),a 



This is necessarily a functional integral, as we must account for a continuum of paths taken by the bj„ (and 
their complex conjugates) as functions of the continuous variable r. Also note, in keeping with convention, we 
refer to the discrete field amplitudes bj a simply as fields. 
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which is clearly proportional to the expansion parameter 1 jz. We exploit this property to 
systematically expand the action up to second order in AS , in order to derive an effective 
action S° s for site 0. The effective action is defined through the corresponding partition 
function Z" ff = Z/Z®, where 

(5) Z e ° ff = JD[b* ]D[b ] exp (-S° ff [b*„ bo]/») , 

(6) Z (0) = rD (0) [b*]D (0) [b]exp(-5 <0) [b*,b]/s). 

In Eq. (|6]l the integral specifically excludes the site fields , as indicated by the notation 
Z) (0) . This is equivalent to performing all integrals appearing in the definition of Z, except 
those located at site 0. In this way we obtain 

z e° ff =i ro[bS]D[bo]exp(-s [b5,bo]/fi) 

(7) J 

x I D (0) [b*]Z) (0) [b]exp(-{5 (0) [b*,b]+A5[b*,b]}/^), 

which can be expanded in powers of AS through 

(8) |D (0) [b*]D (0) [b]exp(-{s (0) [b*,b] + AS[b*,b])/h) 

°° 1 f 
= Z Wk\ J Di0) ^ D(0) W(- AS D»*. b])* exp (-SolK, bol/h) . 

Using now the definition for an expectation value in the cavity system 

(A)(0) = (Z^y 1 f ^>*]D<>]A[b*,bK s< ° ,/; \ 

this immediately gives us 

D[bg]D[bo]exp(-S [b5,bo]/fi)^ , /(0) . 

k=0 

Note that the terms ({AS )*)(0) still depend on the fields at site 0, and cannot be taken outside 
the integral. We determine explicitly the two lowest orders in AS , i.e., to the order in AS 
we wish to consider: 

(10) <AS) ( o) = - - f dr J] J' a [b; )a (T)(b ja (T)} (()) + c.c] , 

Z J ° (0/>,a 

<(A5) 2 ) (0 )=4 [( dTdr'Y Yj'j;[b; )a (T)b; )7 (T')(bj a (T)bj ly (T')} i()) 

(11) Z JJ ° W),aW),y 

+ b* )a (T)bo 7 (T')(b ja (T)by y (T')} {Q) + C.C.]. 

The final (re-exponentiation) step is to note that, taking the effective action to be 

(12) ^ ff [b ,b ] = S [b ,b ] + <AS> ( o) - ^ [<(AS) 2 ) (0) - <AS) 2 0) ] , 

it follows that taking exp(-5° ff [bJ,bo]/^) to second order in AS within Eq. Q yields 
Eq. (|9} to second order in AS (and consequently in 1 /z), as desired. The second-order terms 
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can be conveniently expressed in terms of the connected Green's function in the cavity sys- 
tern, through (AS >* 0) -<(AS ) 2 > (0) = jjf drdr' Z m , a E<o/>, 7 ^^(^(T' T ')B y (r')], 
where B; T (r) = (^(t), b 0a (r)) and 

G >;/r (T ' T) °>;/y (T ' r) 
G (0)«» (t , r) G (0)d ( . , 
j> r j a y ' / frJ<* ' y - 

In the above, we have Gg* y (T, O = <Mt)>(0)<&> 7 (tO>(0)-<M^} 7 (t')>(0), and G^ r (T, r') = 
(bj a (T))(0)(bfy(T'))(0) - (bj a (T)bj>y(T'))(0). Note that the Green's function has a matrix form 
in spin-space as well as Nambu (particle-hole) space; the latter is due to possible off- 
diagonal long-range superfluid order. 

Since the action is invariant under imaginary-time translations single particle expecta- 
tions values like (bj a (T)) do not depend on r. Correlation functions only depend on the 
imaginary -time difference t - t' . This latter fact is used in the following to write the 
equations in terms of the bosonic Matsubara frequencies a>„ = 27m/h/3, reciprocal to this 
imaginary-time difference. 

2.2. Self-Consistency Relations. To obtain a closed self-consistency loop, expectation 
values in the cavity system (i.e., (0)(O)), must be identified with expectation values on the 
impurity site (i.e., (O)o). In the case of the BDMFT equations, this step requires a careful 
treatment, as possible 1 jz corrections must be taken into account. 

We first consider the second-order terms. Because terms to second order in AS within 
the effective action already appear at second order in 1 jz, any 1 /z corrections lead to an 
irrelevant contribution and we need consider leading order behaviour only. We can there- 
fore apply similar reasoning as for the fermionic case [21], and take the limit of infinite 
dimensions in dealing with these terms. 

We combine all terms within the effective action quadratic in the impurity site fields 
into the Weiss field: 

(14) @ a \(iu n ) - 6ay(iha)„cr z + fi a I 2 ) - 2_j "^T" G /"!; y ^ w ")' 

where I2 is the 2x2 identity matrix. From the fact that the self-energy is a local (i.e., 
momentum-independent) quantity in infinite dimensions GTI . then follows the local Dyson 
equation 

(15) G~Jy{i<*>n) = G~ r (/w„) + H ay (ioj n ), 

For future use we also define the hybridization function A oy (ia> n ) = 8 a y(ifko n (T z + Ha^i) ~ 
G~ y (/w„) - 'L a y(iu>„), such that Q ay (ia>„) = 5 ay (itia> n cr z + ii a \2) - A ay (ito n ), and also 

t' J' 

(16) A ay (iaj„) = 2_j ^V" G >;/ y ( /w ")- 

«y>,<o/> Z 

The superfluid order parameter (b a ), on the other hand, appears at first order in 1 jz in 
the effective action. To be consistent, we must therefore take possible 1 jz corrections into 
account, which occur because we must calculate the expectation value of b a in the cavity 
system. We now show that this indeed leads to an important correction. Qualitatively 
this can be understood as being because the sites j on which (b a ) is calculated have one 
fewer neighbour in the cavity system, as the impurity site has been removed. Previously 
this motivated us to implement the 1/z correction by first order perturbation theory in 
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the missing neighbour, which on the Bethe lattice gives rise to results for the superfluid- 
insulator transition very close to the numerically exact solution lfTTll23l . Here we show 
that we can also derive a closed expression without invoking perturbation theory. 

First of all we remark that in the original homogeneous system the expectation value 
is independent of the lattice position and in particular is equal to that at site (chosen 
as the impurity site): (bj a ) = (bo a ). Secondly, expectation values on the impurity site are 
calculated including all relevant orders in 1 jz and thus constitute the DMFT-approximation 
for the expectation values in the original system: <(9o)o = (Go). For the expectation values 
in the cavity system this is different, because the impurity site is missing. 

To calculate the expectation value bj a for a site in the cavity system we couple the 
bosonic fields at site j to a generic source term STja E2- Hence, the partition function 

(17) Z[J ja ,J] a ] = rD[b*]D[b] e - S[b%b]/s+ / dT[ ^ (T) ^' (T)+ ^ (T)/ ^ (r)] 

depends on this source, and we obtain the expectation value from the functional derivative 

(18) <Mr» = t4- MZ[J ja ,J*J). 

J jot 

We now carry out the cavity construction with the sources present in order to derive 
ZgnU'jai ^Tjal- Since we are only interested in calculating the superfluid order parame- 
ter, we only keep terms linear in the sources; to calculate 1/z corrections to the Green's 
function one must retain second order terms as well. 

Defining the shorthand S j — f dr[b* (j)Sia(f) + 3* (j)bj a (j)], we obtain: 

ZP eS [J ja ,J] a l ~ j D[b*]D[bK s ™» +s ' 

(19) ~ [fl[b']/)[b]«- (S » [b ^ l+s '° ,|b * b|)/s 

As previously, we then integrate out the fields in the cavity system and re-exponentiate, 
yielding 

<AS> (0) 



ZeffLfjatJja 



*J= f D[b* ]D[b ]e- s ° [b ° M/n 



1 



Ti 



(20) + _L (( A5 ) 2 )(0)+(5Ao) _<^M + . 



■/ 



£)[| ) *]£)|-|j ] e -S°i r [bJ,bo]/S+(Sj> ( „ ) -«S J AS>,o ) -(5 J >,oi(AS> ( o))//i > 



for which we require the following expectation values: 

-Ti/S 



dT [{bj a (T)r (0) J ja (T) + J* a (T){b ja (T)) (0) ] , 

(SjAS) m =- J] - [[ dTdT'{T ja {T)[bl y (T'){by y (T')b ja (T)\ 
<0/>,y Z JJo ^ 

b 0y (T')(by y (j')bj a (j)\o)\ 



(22) W) ' r 



+ c.c 



M. SNOEK AND W. HOFSTETTER 



Substituting these together with Eq. ( pL0| > into Eq. ( |20l >, we then use Eq. ( [T8| >, and determine 
the following result for the superfiuid order parameter: 



1 / r"p ( 

(b ja (T)) =(b Ja (T)) (0 ) + t V - dr'l 

n <m.y z Jo [ 

( 23 ) x <^ r (r')> [{b ry {T')b ja {T)\ Q) - {b ry {r')) m {b ja {T)) m \ 

+ <V(t')> [(b} y (T')b ja (T)) {0) - <&} r (T')>(0)<MT)>(0)] 

which can be expressed in terms of the Green's function [Eq. ( p"3j )] as: 

(24) (bjjT)} = <^ a (T)> (0 ) 



Using now that the superfiuid order parameter is independent of r, we can carry out the 
integral over r' and obtain 

/' 

(25) (b ja ) = {b ja ) (0) ~ Yj~ [< b U G t°j'y {tD " = 0) + <V><W y (a,„ = 0)] , 

<0/>,y Z 

which indeed yields a 1/z correction to (bj a )(p) with respect to (bj a ). To complete this we 
note that that the source field coupling to K in the effective action is given by J a Ti(0j)(bja)(0)- 
This we can now re-express as 

Ja Z<M(0 ^ J] (bja) + Z 7 [^M'y^ + <V>G^, y (0)] 
<0j) (0» { <0/>,y j 

(26) =^«+ Z ^( G &rW +G 2k<°>) 

<0/>,<0j>,T 



=z/^ a + J] r [a^(0) + a° t (0)] = z/„c, 



where we have introduced </> a - (bj a ) (which is y-independent and assumed to be real) and 
expressed the answer in terms of diagonal and off-diagonal elements of the hybridization 
function at zero frequency. This constitutes the self-consistency equation for the superfiuid 
order parameter. 

For the special case of the Bethe lattice, the hybridization function is proportional to 
the Green's function. In this case, the expression derived from perturbation theory fTT| 
coincides with the closed expression in Eq. ( |26j >, explaining why the exact solution of the 
phase diagram for the single component bosons was so closely reproduced IfTTI . 

2.3. Implementation. In order to complete a self-consistency loop, one has to extract the 
self-energy and the superfiuid order parameter from the effective action such that a new 
input-superfluid order parameter and hybridization function can be calculated. We use 
the exact diagonalization method to do this, which means that the effective action is first 
mapped to a bosonic Anderson Hamiltonian. The Anderson Hamiltonian is not uniquely 
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defined, and we make here the specific choice: 

^And =~ / ; /"««« +A U ay h a [hy - S a y) 
a ay 

(27) 

I \ a J a 

Here the b^ (n a = b a b a ) operators act on the impurity site. In addition we have introduced 
(annihilation) operators af which act on orbitals labelled by I. These orbitals provide a 
non-interacting bath for the impurity site, which incorporates the hybridization with the 
other lattice sites. We have introduced a vector and matrix notation: b a - (b a , b a ) T and 
a/ = (S/, a\) T , and V/ ff = ( w" v '° ) and e; = ( % e j 2 ) are matrices in Nambu space. The 
matrix e\ includes the on-site terms on orbital / (note that there is no dependence on the spin 
index cr\), whereas W contains the diagonal and off-diagonal amplitudes for particles with 
spin index <x to tunnel to orbital /. Integrating out the fields on the orbitals a/, one obtains 
a single-site action with hybridization function: 

(28) A£"W) = I J] Yfa(kW* " e 'y l V 'r 

The parameters in the matrices V/„ and 6/ are obtained by fitting the Anderson hybridiza- 
tion function to the hybridization function obtained from the local Dyson equation. The 
Anderson Hamiltonian can then straightforwardly be implemented in the Fock basis. Af- 
ter diagonalization, the local Green's function can be obtained from the eigenstates and 
eigenenergies: in the Lehmann-representation we obtain 

(29) G a y(ico n ) =- \{m\b a \n){n\by\m)——- — +/3<f> a <f> 7 , 

mn " '" " 

(30) G ay (ito„) =- \(m\b a \n}(n\by\m)——- — +ft a <p r . 

£j Hj II I— 'til "T" LllClJti 

mn 

Knowing the Green's function, the self-energy can in principle be determined from the 
local Dyson equation. However, this can lead to inaccurate results. A more accurate way 
is to use the equations of motion for the Green's function, in the same way as proposed for 
the fermionic case [24|. This leads to the expression: 

(31) l, a y(iu> n ) = 2_j F yay Gy r - zJ a ^ ®ay Gy y 

with <D„ y = (p£<py (\\)<5 n0 and correlation functions F d yay = Ggtg^gt and F° yay = Ggt^gtgt 
[where we note G^^(t) = (T t A(t)B(0)) and T T is the imaginary-time ordering operator], 
which are combined in the matrix F yQ , y according to F = ( ,JL, ( j£., ). 

3. Validity Issues 

3.1. Validity Domain. The BDMFT equations derived here rely on the limit of high di- 
mensions, in particular implying a local self-energy. This means that phases in which the 
self-energy has a non-trivial dependence on momentum (like p-wave or <i-wave pairing) 
cannot be captured. Whereas local correlations are taken into account exactly, non-local 
correlations are treated on the mean-field level. Apart from these restrictions, BDMFT 
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Figure 1 . Phase diagrams obtained by BDMFT for a Bose-Bose mixture 
on the Bethe lattice at total filling n — 1 and z — 4 ifTTIl . 



can be applied in a wide variety of settings; both the homogeneous case and the inhomo- 
geneous case can be treated, the latter by a straightforward extension towards Real-Space 
BDFMT, analogous to Real-Space DMFT for fermionic atoms 12511261 . The scheme is 
completely flexible regarding the spin of the atoms involved [27|. 

3.2. Relevance to Other Theories. The BMDFT equations as derived here incorporate 
Gutzwiller mean-field theory as the lowest order terms in 1/z |frJl[T4l . The second order 
terms correspond to the dynamical mean-field terms familiar from fermionic DMFT [21]. 



4. Application: Ttwo-Component Mixture on the Bethe Lattice 

Applying BDMFT to a Bose-Bose mixture (denoted by b and d) yields a rich phase 
diagram, as shown in Fig. [T] Here we have chosen symmetric interactions (Ui,b = Udd — 
U) and considered the case that Uyi •* U. The total filling is fixed to one: rib + rid = 1. 
The phase diagram includes a superfluid state and insulating states with spin order. The 
spin order is ferromagnetic in the xy-plane when the tunneling ampitudes %j. for the two 
species are comparable and antiferromagnetic in the z -direction for a larger imbalance in 
effective mass [TD] [TT] . The transition between those two types of spin ordered states 
is of first order. For very anisotropic tunneling amplitudes, we find a phase in which the 
species with the small tunneling amplitude becomes localized at every other lattice site, 
whereas the second species is still superfluid IflOl [TT1 [121 . Since in this phase off-diagonal 
superfluid order occurs together with spontaneous translational symmetry breaking, this is 
a supersolid phase IfTTl . 

Above the critical temperature for spin order (cf. the lowest two panels in Fig. fTJ 
the magnetic order disappears, giving rise to an unordered insulator. The supersolid phase 
melts into a monofluid phase, in which one of the species is still superfluid, but translational 
symmetry is restored. 
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